/*
 * C Implementation: ensemble_mean
 *
 * Description: index a GRIB file, select a specific parameter and compute ensemble mean.
 *  索引 GRIB 文件，选择特定参数并计算整体均值
 *
 */
#include "eccodes.h"

int sfx_grib_ensemble_index()
{
    int err = 0;

    char* infile = "d:\\temp\\era5-levels-members.grib";

    // create index of file contents for paramId and number
    // keys 是以逗号分隔的索引键列表 (自定义 && 不能有空格)
    // 后缀表示类型
    // :l for long (or :i)
    // :d for double
    // :s for string
    char* keys         = "paramId:s,number:i";
    codes_index* index = codes_index_new_from_file(0, infile, keys, &err);
    CODES_CHECK(err, 0);

    // ----------------------------------------------------
    // #01: "paramId"
    // get size of "paramId" list
    size_t paramIdSize = 0;
    CODES_CHECK(codes_index_get_size(index, "paramId", &paramIdSize), 0);
    printf("grib contains %zu different parameters\n", paramIdSize);

    // allocate memory for "paramId" list & get list of "paramId"
    char** paramId = (char**)malloc(paramIdSize * sizeof(char*));
    CODES_CHECK(codes_index_get_string(index, "paramId", paramId, &paramIdSize), 0);  // codes_index_get_string (why is string?)
    for (int i = 0; i < paramIdSize; i++) {
        printf("paramId[%d] = %s\n", i, paramId[i]);
    }

    // ----------------------------------------------------
    // #02: "number"
    // get size of ensemble number list
    size_t numberSize = 0;
    CODES_CHECK(codes_index_get_size(index, "number", &numberSize), 0);
    printf("GRIB contains %zu different ensemble members\n", numberSize);

    // allocate memory for ensemble number list & get list of ensemble numbers
    long* number = (long*)malloc(numberSize * sizeof(long));
    CODES_CHECK(codes_index_get_long(index, "number", number, &numberSize), 0);
    for (int i = 0; i < numberSize; i++) {
        printf("number[%d] = %d\n", i, number[i]);
    }

    // ----------------------------------------------------
    // select T850 with paramId 130 (why 130 对应 T850 => 是 grib 规定的?)
    CODES_CHECK(codes_index_select_string(index, "paramId", "130"), 0);

    /** 
     * paramIdSize = 2  : 有2个参数. paramId = ["129", "130"]
     * numberSize = 10  : 每个参数有10组消息 (message subset). number = [0,1,2,3,4,5,6,7,8,9] (数组元素是 message 的索引?)
     * values_len = n   : 每个消息有 n 个数值
     */


    // loop over all members
    int count         = 0;
    double* result    = NULL;
    size_t values_len = 0;
    for (size_t i = 0; i < numberSize; i++) {
        count++;

        // select an individual ensemble number (独立的 number 数据集)
        // 所选子集 (message subset) is into index.
        CODES_CHECK(codes_index_select_long(index, "number", number[i]), 0);

        // create handle for next GRIB message from index
        codes_handle* h = codes_handle_new_from_index(index, &err);
        if (err) {
            fprintf(stderr, "Error: %s\n", codes_get_error_message(err));
            break;
        }

        // get the size of the values array
        CODES_CHECK(codes_get_size(h, "values", &values_len), 0);

        // allocate memory for the GRIB message
        double* values = (double*)malloc(values_len * sizeof(double));

        // allocate memory for result
        if (i == 0) {
            result = (double*)calloc(values_len, sizeof(double));
        }

        // get data values
        CODES_CHECK(codes_get_double_array(h, "values", values, &values_len), 0);

        // add up values
        for (int j = 0; j < values_len; j++)
            result[j] = result[j] + values[j];

        // free memory and GRIB message handle
        free(values);
        CODES_CHECK(codes_handle_delete(h), 0);
    }

    printf("We considered %d ensemble members\n", count);

    for (int j = 0; j < values_len; j++) {
        result[j] = result[j] / count;
    }

    double min = 1e13, max = -1e13, avg = 0;
    for (int j = 0; j < values_len; j++) {
        if (min > result[j]) {
            min = result[j];
        }
        if (max < result[j]) {
            max = result[j];
        }
        avg += result[j];
    }

    avg = avg / values_len;

    printf("==============================================================================\n");
    printf("Stats for ensemble mean of T850\n");
    printf("Min: %f Max: %f Avg: %f\n", min, max, avg);
    printf("==============================================================================\n");

    // free memory
    for (int i = 0; i < paramIdSize; i++)
        free(paramId[i]);
    free(paramId);
    free(number);
    free(result);

    codes_index_delete(index);

    return 0;
}

/*
int main(int argc, char* argv[])
{



    // /* finally free all other memory 


    return 0;
}

*/